function J = Jaco(x, y, z, a, b, c, delta)
    f0 = zeros(3,1);
    [f0(1), f0(2), f0(3)] = SIMM(x, y, z, a, b, c);
    J = zeros(3,3);
    for j = 1:3
        x_perturb = [x; y; z];
        x_perturb(j) = x_perturb(j) + delta;
        [fp1, fp2, fp3] = SIMM(x_perturb(1), x_perturb(2), x_perturb(3), a, b, c);
        J(:, j) = ([fp1; fp2; fp3] - f0) / delta;
    end
end
